# -*- coding: utf-8 -*-
'''
Utilities to plot gps data

@author: Juan C. Espinoza
@contact: jucar.espinoza@gmail.com

December 09, 2009

CHANGELOG:

version 1
- Initial release
'''

__version__ = 1.0

import os, re, sys, time, math
import StringIO
import numpy as np
import Utility as utl
from gpsdatetime import GPSDateTime, datetime, timedelta

import Image
from GifImagePlugin import getheader, getdata

deg2rad = utl.deg2rad
PI      = utl.PI
set     = utl.set
PLT     = False
MPL     = False

def setup(environ=None):
    '''
    '''
    global PLT, MPL
    if environ=='web':
        os.environ['HOME'] = '/var/www/mpl'
        import matplotlib as MPL
        MPL.use('Agg')
    else:
        import matplotlib as MPL
    
    import matplotlib.pyplot as PLT

class MyFigure(object):
    '''
    '''
    def __init__(self, nrows=1, ncols=1, figname=None, **kwargs):
        '''
        '''
        global PLT, MPL

        if not MPL:
            setup(environ=figname)
        
        self.LC = MPL.collections.LineCollection
        self.CM = MPL.cm
        self.py = PLT
        self.figure  = PLT.figure(**kwargs)
        self.figname = figname
        self.axes    = []
        self.nrows   = int(nrows)
        self.ncols   = int(ncols)
        self.rowcnt  = 1
        self.colcnt  = 1

    def adjust(self, left=None, bottom=None, right=None, top=None, wspace=None,
               hspace=None):
        '''
        Subplots adjust in figure
        '''
        width  = self.figure.get_figwidth()
        height = self.figure.get_figheight()
        if not bottom:  bottom  = 0.5*height**-1
        if not top:     top     = 0.74*height**0.11
        if not left:    left    = 0.125
        if not right:   right   = 0.9
        if hspace is None:  hspace  = 0.2
        if wspace is None:  wspace  = 0.2
        self.figure.subplots_adjust(left, bottom, right, top, wspace, hspace)

    def suptitle(self, s, **kwargs):
        '''
        '''
        x = kwargs.pop('x', 0.5)
        y = kwargs.pop('y', 0.98)
        if ('horizontalalignment' not in kwargs) and ('ha' not in kwargs):
            kwargs['horizontalalignment'] = 'center'

        if ('verticalalignment' not in kwargs) and ('va' not in kwargs):
            kwargs['verticalalignment'] = 'top'

        t = self.figure.text(x, y, s, **kwargs)
        return t

    def text(self, x, y, s, **kwargs):
        '''
        '''
        t = self.figure.text(x, y, s, **kwargs)
        return t

    def show(self, watermark=True):
        '''
        '''
        if watermark:
            self.figure.text(1, 0, 'LISN\nLow-Latitude Ionospheric Sensor Network',
                             fontsize=8, color='gray', ha='right', va='bottom',
                             alpha=0.5)
        if self.figname == 'web':
            imgdata = StringIO.StringIO()
            self.figure.savefig(imgdata)
            imgdata.seek(0)
            return imgdata.read()
        elif self.figname:
            self.figure.savefig(self.figname)
            print 'figure saved: %s' % self.figname
        else:
            PLT.show()

    def add_axes(self, shared=False, **kwargs):
        '''
        '''
        polar = kwargs.pop('polar', False)
        axisbg = kwargs.pop('axisbg', None)
        if self.colcnt%(self.ncols+1)==0:
            self.colcnt = 1
            self.rowcnt += 1
        
        ax = MyAxes(**kwargs)
        pos = self.colcnt+((self.ncols)*(self.rowcnt-1))
        ax.n = pos
        ax.ax = self.figure.add_subplot(self.nrows, self.ncols, pos, polar=polar, axisbg=axisbg)
        if shared:
            ax.shax = self.figure.add_axes(ax.ax.get_position(), sharex=ax.ax,
                                           frameon=False)
            ax.shax.yaxis.tick_right()
            ax.shax.yaxis.set_label_position('right')
            ax.ax.yaxis.tick_left()
        self.axes.append(ax)
        self.colcnt += 1
        return ax

    def set_labels(self, **kwargs):
        '''
        '''
        for ax in self.axes:
            ax.set_labels(**kwargs)

    def set_xticks(self, hidden=[], **kwargs):
        '''
        '''
        for ax in self.axes:
            if ax.n in hidden:
                ax.set_xticks(True, **kwargs)
            else:
                ax.set_xticks(**kwargs)

    def set_yticks(self, hide_labels=False, **kwargs):
        '''
        '''
        for ax in self.axes:
            if hide_labels:
                ax.set_yticks(self.ncols, **kwargs)
            else:
                ax.set_yticks(**kwargs)

class MyAxes(object):
    '''
    '''
    def __init__(self, **kwargs):
        self.n            = 1
        self.title        = kwargs.get('title', '')
        self.text         = kwargs.get('text', '')
        self.xlabel       = kwargs.get('xlabel', '')
        self.ylabel       = kwargs.get('ylabel', '')
        self.y2label      = kwargs.get('y2label', '')
        self.xticks       = kwargs.get('xticks', [])
        self.xticklabels  = kwargs.get('xticklabels', [])
        self.yticks       = kwargs.get('yticks', [])
        self.yticklabels  = kwargs.get('yticklabels', [])
        self.y2ticks      = kwargs.get('y2ticks', [])
        self.y2ticklabels = kwargs.get('y2ticklabels', [])
        self.ylim         = kwargs.get('ylim', None)
        self.y2lim        = kwargs.get('y2lim', None)
        self.xlim         = kwargs.get('xlim', None)
        self.xminmax      = set()
        self.yminmax      = set()
        self.y2minmax     = set()
        self.ax           = None
        self.shax         = None

    def set_labels(self,  **kwargs):
        '''
        '''
        if 'size' not in kwargs or 'fontsize' not in kwargs:
            kwargs['fontsize'] = 10
                
        if self.ylabel: self.ax.set_ylabel(self.ylabel, **kwargs)
        if self.xlabel: self.ax.set_xlabel(self.xlabel, **kwargs)
        if self.title : self.ax.set_title(self.title, **kwargs)
        
        if self.shax:                        
            if self.y2label: self.shax.set_ylabel(self.y2label, **kwargs)
            
        if self.text:
            xlim = self.ax.get_xlim()
            ylim = self.ax.get_ylim()
            self.ax.text(xlim[0]+(xlim[1]-xlim[0])*0.95, ylim[1]*0.8,
                         self.text, color='0.4', size=8, ha='right', va='top')

    def set_xticks(self, hidden=False, **kwargs):
        '''
        '''
        
        if self.xticks:
            self.ax.set_xticks(self.xticks)
        if hidden:
            ticklabels = self.ax.get_xticklabels()
            PLT.setp(ticklabels, visible=False)
            self.ax.set_xlabel('')
        else:
            if self.xticklabels:
                self.ax.set_xticklabels(self.xticklabels, **kwargs)
            elif self.xticks:
                self.ax.set_xticklabels(self.xticks, **kwargs)

        if self.xlim  : self.ax.set_xlim(*self.xlim)

        if self.shax:
            ticklabels = self.shax.get_xticklabels()
            PLT.setp(ticklabels, visible=False)

    def set_yticks(self, ncols=None, **kwargs):
        '''
        '''
        if self.yticks:
            self.ax.set_yticks(self.yticks)

        if self.yticklabels:
            self.ax.set_yticklabels(self.yticklabels, **kwargs)
        elif self.yticks:
            self.ax.set_yticklabels(self.yticks, **kwargs)

        if self.ylim  : self.ax.set_ylim(*self.ylim)
        
        if self.shax:
            if self.y2ticks:
                self.shax.set_yticks(self.y2ticks)
            if self.y2ticklabels:
                self.shax.set_yticklabels(self.y2ticklabels, **kwargs)
            elif self.y2ticks:
                self.shax.set_yticklabels(self.y2ticks, **kwargs)

            if self.y2lim  : self.shax.set_ylim(*self.y2lim)

        if ncols:
            cols = range(ncols)
            if self.col in cols[:-1] and self.shax:
                self.shax.set_ylabel('')
                ticklabels = self.shax.get_yticklabels()
                PLT.setp(ticklabels, visible=False)
            if self.col in cols[1:]:
                self.ax.set_ylabel('')
                ticklabels = self.ax.get_yticklabels()
                PLT.setp(ticklabels, visible=False)

    def plot(self, *args, **kwargs):
        '''
        '''

        shared = kwargs.pop('shared', False)

        if shared:
            return self.shax.plot(*args, **kwargs)
        else:
            return self.ax.plot(*args, **kwargs)

class myFigure(object):
    '''
    '''
    def __init__(self, figname=None, **kwargs):
    
        err1 = '''Please create a full access folder for Matplotlib 
        Environment vars for example /var/www/mpl/'''
        
        if figname == 'web':
            os.environ['HOME'] = '/var/www/mpl'
        import matplotlib
        if figname:
            matplotlib.use('Agg')
        
        import pylab
        from matplotlib.collections import LineCollection
        self.py = pylab
        self.LC = LineCollection
        self.figure = pylab.figure(**kwargs)
        self.figname = figname
        
    def adjust(self, left=None, bottom=None, right=None, top=None, wspace=None,
               hspace=None):
        '''
        Subplots adjust in figure
        '''
        width  = self.figure.get_figwidth()
        height = self.figure.get_figheight()
        if not bottom:  bottom  = 0.5*height**-1
        if not top:     top     = 0.74*height**0.11
        if not left:    left    = 0.125
        if not right:   right   = 0.9
        if not hspace:  hspace  = 0.2
        if not wspace:  wspace  = 0.2
        self.figure.subplots_adjust(left, bottom, right, top, wspace, hspace)
        
    def show(self):
        '''
        '''
        self.figure.text(1, 0, 'LISN\nLow Latitude Ionospheric Sensor Network', 
                         fontsize=8, color='gray', ha='right', va='bottom',
                         alpha=0.5)
        if self.figname == 'web':
            f = os.tmpfile()
            self.figure.savefig(f)
            f.seek(0)
            rawfig = f.read()
            f.close()
            return rawfig
        elif self.figname:
            self.figure.savefig(self.figname)
            print 'figure saved: %s' % self.figname
        else:
            self.py.show()
            
class images2Gif(object):
    '''
    Provides a function (writeGif) to write animated gif from a series
    of PIL images or numpy arrays.

    This code is provided as is, and is free to use for all.

    Almar Klein (June 2009)

    - based on gifmaker (in the scripts folder of the source distribution of PIL)
    - based on gif file structure as provided by wikipedia
    '''
    
    def __init__(self, fname, images, durations=None, loops=0):
        """Given a set of images writes the bytes to the specified stream. """
        fp = open(fname, 'wb')
        if not durations:
            durations = [0.5 for im in images]
        frames = 0
        previous = None
        for im in images:
            if not previous:
                palette = getheader(im)[1]
                data = getdata(im)
                imdes, data = data[0], data[1:]
                header = self.getheaderAnim(im)
                appext = self.getAppExt(loops)
                graphext = self.getGraphicsControlExt(durations[0])
                # write global header
                fp.write(header)
                fp.write(palette)
                fp.write(appext)
                # write image
                fp.write(graphext)
                fp.write(imdes)
                for d in data:
                    fp.write(d)
            else:
                # gather info (compress difference)              
                data = getdata(im)
                imdes, data = data[0], data[1:]
                graphext = self.getGraphicsControlExt(durations[frames])
                # write image
                fp.write(graphext)
                fp.write(imdes)
                for d in data:
                    fp.write(d)
            previous = im.copy()
            frames = frames + 1
        fp.write(";")  # end gif
        print frames, 'frames written'
        fp.close()

    def intToBin(self, i):
        """ Integer to two bytes """
        # devide in two parts (bytes)
        i1 = i % 256
        i2 = int( i/256)
        # make string (little endian)
        return chr(i1) + chr(i2)

    def getheaderAnim(self, im):
        """ Animation header. To replace the getheader()[0] """
        bb = "GIF89a"
        bb += self.intToBin(im.size[0])
        bb += self.intToBin(im.size[1])
        bb += "\x87\x00\x00"
        return bb


    def getAppExt(self, loops=0):
        """ Application extention. Part that secifies amount of loops. 
        if loops is 0, if goes on infinitely.
        """
        bb = "\x21\xFF\x0B"  # application extension
        bb += "NETSCAPE2.0"
        bb += "\x03\x01"
        if loops == 0:
            loops = 2**16-1
        bb += self.intToBin(loops)
        bb += '\x00'  # end
        return bb

    def getGraphicsControlExt(self, duration=0.1):
        """ Graphics Control Extension. A sort of header at the start of
        each image. Specifies transparancy and duration. """
        bb = '\x21\xF9\x04'
        bb += '\x08'  # no transparancy
        bb += self.intToBin( int(duration*100) ) # in 100th of seconds
        bb += '\x00'  # no transparant color
        bb += '\x00'  # end
        return bb

def png2gif(fname, pngs, resize=None):
    '''
    '''
    print 'creating GIF...' 
    images = [Image.open(png) for png in pngs]
    if resize:
        images = [img.resize(resize, Image.ANTIALIAS) for img in images]
    images = [im.convert('P',dither=1) for im in images]
    convert = images2Gif(fname, images)

def plot_tec(gdo, figname=None, localtime=False, prns=None, min_ele=25):
    '''
    '''
    
    if localtime:
        timezone = timedelta(gdo.tzoffset/24.)
    prns = [x for x in gdo.prns if x in gdo.phasearcs and x not in gdo.bad_prn]
    prns.sort()
    fig = MyFigure(figname=figname, figsize=(8,6))
    ax = fig.figure.add_subplot(1,1,1)
    p_color = [[ob[0] for ob in gdo.iterlist(prn, ('lat'), ('ele', min_ele))] for prn in prns]
    pmax = max([max(x) for x in p_color if len(x)<>0])
    pmin = min([min(x) for x in p_color if len(x)<>0])
    maxtec = [[ob[0] for ob in gdo.iterlist(prn, ('eqTEC'), ('ele', min_ele))] for prn in prns]
    maxtec = max([max(x) for x in maxtec if len(x)<>0])
    mintec = [[ob[0] for ob in gdo.iterlist(prn, ('eqTEC'), ('ele', min_ele))] for prn in prns]
    mintec = min([min(x) for x in mintec if len(x)<>0])
    dt1 = gdo[-1].epoch
    dt0 = dt1 - timedelta(1)
    
    for prn in prns:
        points = [(dt.doys, ob) for dt, ob in \
                  gdo.iterlist(prn, ('epoch', 'eqTEC'), ('ele', min_ele)) if dt>=dt0]
        if len(points)==0:
            continue
        #ax.annotate(prn, points[0], size=8)
        ax.text(points[0][0], points[0][1], prn, size=8)
        segments = []
        for i in range(len(points)-1):
            if points[i+1][0]-points[i][0] > 1./12:
                #ax.annotate(prn, points[i+1], size=8)
                ax.text(points[i+1][0], points[i+1][1], prn, size=8)
                continue
            else:
                segments.append((points[i],points[i+1]))
        LC = fig.LC(segments, linestyle='dotted')
        LC.set_linewidth(1)
        p_color = [ob[0] for ob in gdo.iterlist(prn, ('lat'), ('ele', min_ele))]
        LC.set_array(np.array(p_color+[pmin, pmax]))
        ax.add_collection(LC)

    cbar = fig.figure.colorbar(LC,orientation='vertical', fraction=0.05,
                               shrink=0.5)
    cbar.set_label('Lat.')#, size=10)
    #cbar.ax.set_xlim(pmin, pmax)
    for label in cbar.ax.get_yticklabels():
        label.set_fontsize(8)

    ax.set_xlim(dt0.doys, dt1.doys)
    ax.set_ylim(0, (int(maxtec)+5)/5*5)
    dt0 = GPSDateTime(dt0.year, dt0.month, dt0.day, dt0.hour)
    dt1 = GPSDateTime(dt1.year, dt1.month, dt1.day, dt1.hour)
    tk_list = [dt0+timedelta(x/24.) for x in range(1,25)]
    ax.set_xticks([dt.doys for dt in tk_list])
    if localtime:
        ax.set_xticklabels([`(dt+timezone).hour` for dt in tk_list])
        ax.set_xlabel('Time [LT]')
    else:
        ax.set_xticklabels([`dt.hour` for dt in tk_list])
        ax.set_xlabel('Time [UT]')
    ax.set_ylabel('TEC [Units]')

    
    if 'fullname' in gdo.station:
        ax.set_title(gdo.station['fullname'], size=12)
    else:
        ax.set_title('Station: %s' % gdo.header['marker'][0].upper(), size=12)
    
    if gdo[0].epoch.date==gdo[-1].epoch.date:
        fig.figure.text(0.05, 0.95, dt1.strftime('%Y/%m/%d'),
                        size=10, ha='left')
    else:
        fig.figure.text(0.05, 0.95, '%s - %s' % (dt0.strftime('%Y/%m/%d'),
                        dt1.strftime('%Y/%m/%d')), size=10, ha='left')
    
    if gdo.rec_bias<>0:
        fig.figure.text(0.95, 0.95, 'rx_bias = %.1f' % gdo.rec_bias,
                        size=10, ha='right')

    return fig.show()

def plot_gpsdata(gdo, figname=None, localtime=False):
    '''
    '''
    if localtime:
        timezone = timedelta(gdo.tzoffset/24.)
    prns = [x for x in gdo.prns if x in gdo.phasearcs]
    prns.sort()
    fig = MyFigure(figname=figname, figsize=(7,7))
    ax = fig.figure.add_subplot(1,1,1)

    dt1 = gdo[-1].epoch
    dt0 = dt1 - timedelta(1)
    data = {}
    all = []
    for rec in gdo:
        time = rec.epoch.doys        
        for prn in prns:
            if prn not in data:
                data[prn] = []
            if prn in rec:
                data[prn].append(time)
        all.append(time)

    prns = data.keys()
    prns.sort()
    i = 2
    for prn in prns:
        ax.plot(data[prn], [i for x in data[prn]], '|', ms=4, color='blue')
        i+=1

    ax.plot(all, [1 for x in all], '|', ms=3, color='green')
    ax.grid(True)

    ax.set_xlim(dt0.doys, dt1.doys)
    ax.set_ylim(0, len(prns))
    ax.set_yticks(range(i+1))
    ax.set_yticklabels(['', 'All']+prns)
    dt0 = GPSDateTime(dt0.year, dt0.month, dt0.day, dt0.hour)
    dt1 = GPSDateTime(dt1.year, dt1.month, dt1.day, dt1.hour)
    tk_list = [dt0+timedelta(x/24.) for x in range(1,25)]
    ax.set_xticks([dt.doys for dt in tk_list])
    if localtime:
        ax.set_xticklabels([`(dt+timezone).hour` for dt in tk_list])
        ax.set_xlabel('Time [LT]')
    else:
        ax.set_xticklabels([`dt.hour` for dt in tk_list])
        ax.set_xlabel('Time [UT] - %s' % gdo.date.strftime('%Y/%m/%d'))
    ax.set_ylabel('PRN')
    ax.set_title('GPS Data', size=12)    
    fig.adjust()
    return fig.show()

def plot_mag(mdo, figname=None, station=None, figure=None, scale=None):
    '''
    '''
    if not figure:
        fig = MyFigure(figname=figname, figsize=(6,3))
    else:
        fig = figure
    
    H = [x-mdo[0]['H'] for x in mdo.iterlist(['H'])]
    Z = [x-mdo[0]['Z'] for x in mdo.iterlist(['Z'])]
    D = [(x-mdo[0]['D'])*180/3.1415 for x in mdo.iterlist('D')]
    T = [x.hours for x in mdo.iterlist(['epoch'])]
    
    if len(H)>0:
        max_value = max(abs(max(H)),abs(min(H)))
    else:
        max_value = 0
    
    if not scale: scale = (int(max_value*10/10))+(int(0.2*(max_value*10/10)))
    else: scale = int(scale)

    [y1max, y1min, y2max, y2min] = [scale, -scale, scale/10., -scale/10.]
    
    ax  = fig.add_axes(shared=True,
                       xlim=(0,24), ylim=(y1min,y1max), y2lim=(y2min,y2max),                       
                       xlabel='Hour [UTC]',
                       ylabel='H & Z [nT]',
                       y2label='D [min]',
                       xticks=range(0,24,2),
                       yticks=[y1min+(y1max-y1min)*x/4 for x in range(5)],
                       y2ticks=[y2min+(y2max-y2min)*x/4 for x in range(5)],
                       )
    msg = '\nNo Data Available'
    if len(mdo)>0:
        for arc in mdo.arcs:
            i,j = arc
            ax.plot(T[i:j], H[i:j], 'b')
            ax.plot(T[i:j], Z[i:j], 'g')
            ax.plot(T[i:j], D[i:j], 'r', shared=True)    
        ax.ax.text(23.5, y1max-0.2*y1max, mdo.date.strftime('%d/%m/%Y'), size=8, ha='right', va='top')
        msg = ''
    
    fig.set_xticks(size=8)
    fig.set_yticks(size=8)
    fig.set_labels()
    fig.text(0.4,0.96,'H---', size=10, color = 'b', ha='center', va='bottom')
    fig.text(0.5,0.96,'Z---', size=10, color = 'g', ha='center', va='bottom')
    fig.text(0.6,0.96,'D---', size=10, color = 'r', ha='center', va='bottom')
    
    if station:
        ax.ax.text(0.5, y1max-0.2*y1max, station['fullname']+msg, size=8, ha='left', va='top')
    ax.ax.grid(color='0.8')
    fig.adjust(right=0.85)

    if not figure:
        return fig.show()
    else:
        return ax

def plot_s4(sdo, plot_type='0', figname=None, station={}, prns=None, min_ele=25,localtime=False):
    '''
    '''
    if not station:
        station['fullname']='Unknown'
    if not prns:
        prns = [prn for prn in sdo.prns]
        prns.sort()

    if localtime and station:
        timezone = timedelta((station['longitude']/15)/24.)
    else:
        localtime = False

    dt1 = sdo[-1].epoch
    dt0 = sdo[0].epoch
    if dt1-dt0<20*60*60 or dt1-dt0>24*60*60:
        dt0 = dt1 - timedelta(1) 

    dt0 = GPSDateTime(dt0.year, dt0.month, dt0.day, dt0.hour)
    dt1 = GPSDateTime(dt1.year, dt1.month, dt1.day, dt1.hour)
    tk_list = [dt0+timedelta(x/24.) for x in range(0,25)]

    if localtime:
        if (dt1+timezone).doy==(dt0+timezone).doy:
            date = (dt0+timezone).strftime('%Y/%m/%d')
        else:
            date = '%s-%s' % ((dt0+timezone).strftime('%Y/%m/%d'), (dt1+timezone).strftime('%Y/%m/%d'))
        xticklabels = [`(dt+timezone).hour` for dt in tk_list]
        xlabel= 'Time [LT]\n%s' % date
    else:
        if dt1.doy==dt0.doy:
            date = dt0.strftime('%Y/%m/%d')
        else:
            date = '%s-%s' % (dt0.strftime('%Y/%m/%d'), dt1.strftime('%Y/%m/%d'))
        xticklabels = [`dt.hour` for dt in tk_list]
        xlabel = 'Time [UT]\n%s' % date


    if plot_type == '0':
        plots = {}      
        fig  = MyFigure(figname=figname, figsize=(6, 3))
        cmap = fig.CM.get_cmap("jet", 32)
        ax  = fig.add_axes(xlim=((sdo[-1].epoch-timedelta(1)).doys, sdo[-1].epoch.doys), 
                           ylim=(0,1),
                           xlabel=xlabel,
                           ylabel='Scintillation index S4',
                           xticks=[dt.doys for dt in tk_list],
                           xticklabels = xticklabels,
                           yticks=[0.2, 0.4, 0.6, 0.8, 1.0],
                           )
        for prn in prns:
            if prn>32: continue
            for arc in sdo.arcs[prn]:           
                tups = [(rec.epoch.doys, rec[prn]['s4']) for rec in sdo[arc[0]:arc[1]] \
                        if rec[prn]['el']>=min_ele]
                x   = [tup[0] for tup in tups]
                y   = [tup[1] for tup in tups]
                p = ax.plot(x, y, color=cmap(prn-1), label='PRN %02d' % prn)
                dum = [val>0.8 for val in y]
                if True in dum and prn not in plots:
                    plots[prn] = p[0]
        
        ax.ax.grid(True)    
        ax.ax.xaxis.grid(False)
        ax.ax.text(tk_list[-1].doys,0.85, r'$S4\ values\ for\ elevation\ angles\ >\ 25\degree$', 
                   ha='right', size='10', color='0.5')
        fig.set_xticks(size=8)
        fig.set_yticks(size=8)
        fig.set_labels()
        fig.adjust(bottom=0.2,left=0.1, right=0.84)
        fig.suptitle('\n%s' % (station['fullname'], ))
        if plots:
            labels = ['PRN %02d' % s for s in plots.keys()]
            leg=fig.figure.legend(plots.values(), labels, loc=7, prop={'size':8}, borderpad=0.5)
            leg.set_title("S4>0.3")

    elif plot_type == '1':
        ###create fig and axes
        fig = MyFigure(figname=figname, figsize=(5, 5.5), facecolor='white')
        ax  = fig.add_axes(xlim=(-100,100), ylim=(-100,100))
        ax.ax.axison = False
        ###plot custom grid
        ax.plot([-100, 100], [0, 0], color='0.6', lw=0.5)
        ax.plot([0, 0], [-100, 100], color='0.6', lw=0.5)
        for loc in [30,60,0]:
            th = np.linspace(0, 2*utl.PI, 100)
            xr = 100*np.cos(loc*utl.deg2rad)*np.cos(th)
            yr = 100*np.cos(loc*utl.deg2rad)*np.sin(th)
            ax.plot(xr,yr,color='0.6', lw=0.5)
        ###plot coloured line for time arc for prn in data
        for prn in prns:
            if prn>32: continue
            for arc in sdo.arcs[prn]:
                s4_array = np.array([rec[prn]['s4'] \
                                    for rec in sdo[arc[0]:arc[1]]]+[0,1])
                el_array = np.array([100*np.cos(rec[prn]['el']*utl.deg2rad) \
                                    for rec in sdo[arc[0]:arc[1]]])
                az_array = np.array([utl.PI/2-rec[prn]['az']*utl.deg2rad \
                                    for rec in sdo[arc[0]:arc[1]]])
                X = el_array*np.cos(az_array)
                Y = el_array*np.sin(az_array)
                points = zip(X,Y)
                segments = zip(points[:-1], points[1:])
                LC = fig.LC(segments)
                LC.set_linewidth(2)
                LC.set_array(s4_array)
                ax.ax.add_collection(LC)
                ax.ax.text(X[-1],Y[-1], prn, ha='center', va='center',
                           fontsize=9, color='0.5')
        ###add colorbar and labels to figure
        cbar = fig.figure.colorbar(LC, orientation='horizontal', ax=ax.ax,
                                   shrink=0.8, pad=0.08, fraction=0.05)
        cbar.set_label('GPS S4', size=10)
        ticklabels = cbar.ax.get_xticklabels()
        fig.py.setp(ticklabels, size=10)
        ax.set_labels()
        ax.ax.text(102, 0,'E', va='center')
        ax.ax.text(-102,0,'W', va='center', ha='right')
        ax.ax.text(0, 102,'N', ha='center')
        ax.ax.text(0,-102,'S', ha='center', va='top')
        fig.suptitle('GPS S4 in Skymap at %s\n%s ' % (station['fullname'],
                                                      sdo.date))

    return fig.show()

def plot_lines(data, fmt=[], figname=None, title='', **kwargs):
    '''
    Plot x vs a, x vs b etc. in one figure where:
    
    x: datetime serie
    
    data = [(x1, a1, b1...),(x2, a2, b2...),...]
    fmt  = [(var_position, line_format, y_label, comment),...]
    
    '''
    
    date = data[0][0].date()
    x_list = [bub[0].hours for bub in data]
    
    if not fmt:
        c = range(1,len(bubs[0]))
        m = ['-' for x in range(len(fmt))]
        l = ['' for x in range(len(fmt))]
        data_info = zip(c,m,l,l)
    
    kwargs['figsize'] = (6, 0.8*len(fmt)+3.2)
    
    fig = myFigure(figname, **kwargs)
    fig.adjust(hspace = 0.1)

    for i, tup in enumerate(fmt):
        ax = fig.figure.add_subplot(len(fmt), 1, i+1)
        Y = [bub[tup[0]] for bub in data]
        ax.plot(x_list, Y, tup[1])
        yticklabels = ax.get_yticklabels()
        fig.py.setp(yticklabels, size=10)
        ax.set_xticks(range(0,24))
        ax.set_xlim(0, 24)
        ax.set_ylim(*tup[2])
        ax.set_ylabel(tup[3])
        ax.text(23, 0.7*tup[2][1], tup[4], size=8, ha='right')
        if not (i==len(fmt)-1):
            xticklabels = ax.get_xticklabels()
            fig.py.setp(xticklabels, visible=False)
    xticklabels = ax.get_xticklabels()
    fig.py.setp(xticklabels, size=10)
    ax.set_xticks(range(0,24))
    ax.set_xlabel('Time [UT]')
    fig.figure.suptitle(title + '\n' + \
                 datetime.strftime(date, '%m/%d/%Y'), ha='center', va='top')
    fig.show()

def plot_lines2(data, fmt=[], figname=None, title='', rowcol=(1,1), adjargs={},
                **kwargs):
    '''
    Plot x vs a, x vs b etc. in one figure where:
    
    x: datetime serie
    
    data = [[(x1, a1), (x2, a2),...], [(y1, b1), (y2, b2),...],...]
    fmt  = [(list_position, line_format, y_label, comment),...]
    
    '''

    if not fmt:
        c = range(len(data))
        m = ['-' for x in range(len(data))]
        l = ['' for x in range(len(data))]
        fmt = zip(c,m,l,l)
    if 'figsize' not in kwargs:
        kwargs['figsize'] = (6, 0.8*len(data)+3.2)

    fig = myFigure(figname, **kwargs)
    fig.adjust(**adjargs)
    axes = []
    dt1 = datetime(1900,1,1)
    for i, vars in enumerate(data):
        fig_pos = rowcol+(i+1,)
        ax = fig.figure.add_subplot(*fig_pos)
        axes.append(ax)
        dts = [tup[0] for tup in vars]
        if max(dts)>dt1: dt1 = max(dts)
        X = [x.doys for x in dts]
        Y = [tup[1] for tup in vars]
        ax.plot(X, Y, fmt[i][1], ms=1)
        yticklabels = ax.get_yticklabels()
        fig.py.setp(yticklabels, size=8)
        ax.set_ylim(*fmt[i][2])
        ax.set_ylabel(fmt[i][3])
        xticklabels = ax.get_xticklabels()
        fig.py.setp(xticklabels, visible=False)
    dt1_o = dt1
    dt1 = datetime(dt1.year, dt1.month, dt1.day, dt1.hour)
    dt0 = dt1-timedelta(1)
    max_dt = dt1_o.doys
    tk_list = [dt0+timedelta(x/12.) for x in range(1,13)]
    for i, ax in enumerate(axes):
        ax.set_xlim(max_dt-1, max_dt)
        ax.set_xticks([dt.doys for dt in tk_list])
        ax.set_xticklabels([`dt.hour` for dt in tk_list])
        ax.text(max_dt-0.18, 0.7*fmt[i][2][1], fmt[i][4], size=8, ha='left',
                color='g')
    
    last_axes = [axes[-1*x-1] for x in range(rowcol[1])]
    for ax in last_axes:
        xticklabels = ax.get_xticklabels()
        fig.py.setp(xticklabels, visible=True, size=8)
        ax.set_xlabel('Time [UT]')
    fig.figure.suptitle(title + '\n' + \
               dt1_o.strftime('%Y/%m/%d - %H:%M UT'), ha='center', va='top')
    fig.show()

def plot_lines3(data, figname=None, title='', labels = [], rowscols=(1,1), 
                shareax=[], figargs={}, adjargs={}, axsargs=[], pltargs=[]):
    '''
    data = [[(x1, a1), (x2, a2),...], [(y1, b1), (y2, b2),...],...]
    x      : datetimes objects
    a,b    : integers or floats

    figname = filename for save the figure, if is None the figure is shown.
    title   = title for figure.
    rowscols= numbers of rows and cols for multiple plots.
    shareax = List of True or False to indicate if Y values share x axes.
    ###args = dictionary or list of dictionarys with extra arguments for figure,
              adjust_figure, axes and plots.
    
    '''
    N = len(data)
    if not labels:
        labels = ['' for i in range(N)]
    if not shareax:
        shareax = [False for i in range(len(data[0][0])-1)]
    if not axsargs:
        axsargs = [{} for i in range(len(data[0][0])-1)]
    if not pltargs:
        pltargs = [{} for i in range(len(data[0][0])-1)]
    
    fig = myFigure(figname, **figargs)
    fig.adjust(**adjargs)
    axes_left  = []
    axes_right = []
    x_max = GPSDateTime(1900,1,1)
    for i in range(N):
        ax_pos = rowscols + (i+1,)
        dts = [tup[0] for tup in data[i]]
        if max(dts)>x_max:
            x_max = max(dts)
        X = [x.doys for x in dts]
        ax = fig.figure.add_subplot(*ax_pos, **axsargs[0])
        axes_left.append(ax)
        for j in range(len(data[i][0])-1):
            Y = [tup[j+1] for tup in data[i]]
            if not shareax[j]:
                ax.plot(X, Y, **pltargs[0])
                yticklabels = ax.get_yticklabels()
                fig.py.setp(yticklabels, size=8)
                xticklabels = ax.get_xticklabels()
                fig.py.setp(xticklabels, visible=False)
                if 'ylim' in axsargs[0]:
                    ax.set_ylim(axsargs[0]['ylim'])
            else:
                ax1 = fig.figure.add_axes(ax.get_position(True), sharex=ax,
                                          frameon=False, **axsargs[1])
                ax1.yaxis.tick_right()
                ax1.yaxis.set_label_position('right')
                ax.yaxis.tick_left()
                ax1.xaxis.set_visible(False)
                axes_right.append(ax1)
                ax1.plot(X, Y, **pltargs[1])
                yticklabels = ax1.get_yticklabels()
                fig.py.setp(yticklabels, size=8)
                xticklabels = ax1.get_xticklabels()
                fig.py.setp(xticklabels, visible=False)
                if 'ylim' in axsargs[1]:
                    ax1.set_ylim(axsargs[1]['ylim'])

    dt_tick = GPSDateTime(x_max.year, x_max.month, x_max.day, x_max.hour)
    tk_list = [dt_tick-timedelta(hours=2*x) for x in range(11,-1,-1)]

    for i, ax in enumerate(axes_left):
        ax.set_xlim(x_max.doys-1, x_max.doys)
        ax.set_xticks([dt.doys for dt in tk_list])
        ax.set_xticklabels([`dt.hour` for dt in tk_list])
        ax.text(ax.get_xlim()[1]-0.2, 0.7*ax.get_ylim()[1], 
                labels[i], size=8, ha='left', color='g')

    ## Set visible ticks for last axes
    for ax in axes_left[-1*rowscols[1]:]:
        xticklabels = ax.get_xticklabels()
        fig.py.setp(xticklabels, visible=True, size=8)
        ax.set_xlabel(axsargs[0]['xlabel'])

    fig.figure.suptitle(title+'\n'+x_max.strftime('%Y/%m/%d - %H:%M UT'),
                        ha='center', va='top')

    return fig
    

def plot_polar(data, figname=None, title='', adjargs={}, **kwargs):
    '''
    '''
    fig = myFigure(figname, **kwargs)
    fig.adjust(**adjargs)
    ax = fig.figure.add_subplot(111, polar=True)#, axisbg='#E2EDF9')
    prns = data.keys()
    C = []
    R = []
    T = []
    for prn in prns:
        r = 1-data[prn]['el']/90.
        R.append(r)
        theta = -deg2rad*data[prn]['az']+PI/2
        T.append(theta)
        C.append(data[prn]['s4'])
    R.append(0.95)
    T.append(0)
    C.append(0)
    area = [30*x for x in C]
    sct = ax.scatter(T, R, c=C, s=area)
    sct.set_alpha(0.75)
    for i,prn in enumerate(prns):
        ax.text(T[i], R[i], '%i' % (prn, ), size=8, color='g')
    ax.set_thetagrids([0, 45, 90, 135, 180, 225, 270, 315], ['E', '', 'N', '', 'W', '', 'S', ''])
    ax.set_rgrids([0.0000001, 0.3333333, 0.666666666, 1], ['','','',''])
    fig.figure.suptitle(title + '\n' + \
               datetime.strftime(data.epoch, '%Y/%m/%d - %H:%M UT'), ha='center', va='top', size=10)
    fig.show()
